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Abstract 

An analytical approach is developed for bosonic probability amplitudes in unitary linear net- 
works, such as the optical multiport devices for photons. The approach applies for a large number 
of bosons A^ ^ M in the M-mode network. The probability amplitudes of A^ bosons unitarily 
transformed from the input modes to the output modes of a unitary linear network are expressed 

Q-c through a multidimensional integral with the integrand containing a large parameter (A^) in the 

-(— > I 

Ch . exponent. The integral representation allows an asymptotic estimate for the bosonic matrix per- 
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manents up to a small multiplicative error of order 1/A^. The estimate depends on the solution of 
the scaling problem of the unitary M x M-dimensional network matrix: to find the left and right 
K^ . diagonal matrices which scale the unitary matrix to a matrix which has specified rows and columns 

in 

t^^ I sums (equal to the distributions of bosons in the input and output modes). Such scaled matrices 

^^ ' give the saddle points of the integral. In the case of simple saddle points an explicit formula giving 

^D ' an asymptotic estimate for the bosonic probability amplitudes is given. Good accuracy is shown 

cn 

by testing the approximation in the simplest case of just two- mode network (the beam-splitter), 
where the saddle-points are the roots of a quadratic. 
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I. INTRODUCTION 

Linear optical networks have always been playing an important role for the quantum 
manifestations of light. Indeed, The well-known Hong-Ou-Mandel (HOM) dip [l| (see also 
:lefs. |2|,|3|) is a direct manifestation of the quantum indistinguishability of photons. Recently 
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a generalization of the HOM effect was demonstrated in 



Jae general setting of the 
). There are important 



Bell multiport beam splitters (see also, for instance, Refs. p, |6| 

results on the quantum interference phenomena in the multiport devices for the bosons and 

fermions. For instance, the use of symmetry has allowed to show the existence of the zero 

the boson and fermion statistics in the unitary linear networks were recently discovered. [8[ 
There are important experimental advances in the quantum interference experiments with 
indistinguishable particles in linear multiport devices (see, for instance, the very recent 
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multiphoton interference experiment on an integrated optical device [9|). Moreover, the 
experimental efforts are now also directed at building the so-called boson sampler network 
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- |l2| (see also below). Finally, the linear bosonic network lies at heart of an ingenious 
proposal of quantum computation based on the linear optics [l3|. 

It is known that the outcome probability amplitudes in linear bosonic networks (e.g. 
in the optical multiports) are expressed through the matrix permanents, somewhat similar 
to the fermonic amplitudes which are given by the matrix (a.k.a. Slater) determinants. 
Considering A^ non-interacting bosons in an unitary network of M input and output modes, 
one has to compute the matrix permanents of the N x A^- dimensional complex matrices 
composed of repeated rows and columns of the network unitary matrix to find the bosonic 
transition amplitudes of the network. The permanent of a matrix [ij] can be obtained by 
taking the well-known Laplace expansion formula for the matrix determinant and setting 
to "+" the signatures of all permutations in the Laplace summation. The relation of the 
bosonic amplitudes to matrix permanents was explored a long time ago 15| in connection 
with the quantum fields theory. 

Recently the relation of unitary linear bosonic networks to matrix permanents has at- 
tracted a lot of renewed attention. One of the reasons is that the similarity between fermions 
and bosons does not go along as far as one tries to compute the matrix permanents: unlike for 
the matrix determinants, computation of permanent of an arbitrary matrix is #P-complete, 



which is the result of a classic paper on the complexity theory [16]. This means that no 
algorithm polynomial in the matrix size can compute the matrix permanents. The fastest 
known algorithm for computing the permanent of an arbitrary (complex) A^ x A^- dimensional 
matrix is based on Ryser's formula [17] and requires 0{2^^^N'^) flop operations. An interest- 



ing "physical" proof of the matrix permanent comp 



1, H 



computing proposal, [13| was recently discovered. 



exity result, based on the linear optical 
18( 1 Moreover, even the problem of ap- 



proximating the permanent of an arbitrary complex matrix to a polynomial multiplicative 
error was also shown [l8| to be a ^^P-hard problem.^ 

The complexity of computing the matrix permanent was analyzed also in the context of 



the quantum mechanics for the first time in Ref. [20], where a method based on the quantum 
measurement was proposed to directly measure the matrix permanent in a quantum bosonic 
multiport device. It was shown that the permanent of an arbitrary matrix can be expressed 
as a quantum observable. However, the catch of the method lies in the exponential number 
of the necessary measurements, since the variance of such an observable is exponential in 
the matrix size 2y] (for comparison: the matrix determinant can be found in just a single 
measurement). 

The deep connection between the complexity of bosonic networks and that of the matrix 
permanents was used in the recent proposal of a new model of quantum computer with 
noninteracting identical bosons, which, though not being an universal quantum computer, 
nevertheless can perform the computations considered to be hard on a classical computer. 
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21| Such quantum computer was compared to classical Galton's board where, instead of 
the classical balls and a single entry point, the identical bosons are launched into different 
modes of a linear network. The very crucial difference, however, is that in the quantum 
network case the bosonic probabilities themselves are not known beforehand, since they 
cannot be effectively computed on a classical computer (for a sufficiently large network). One 
has to run the actual sampling experiments to find them. In the technical part the proposal 
depends on the hardness to approximate the matrix permanent of a large A^ x A^- dimensional 
submatrix of arbitrary unitary M x M-dimensional matrix and some numerically tested 
conjectures. This proposal generated the experimental efforts on building the necessary 



^ There is an important exception of the permanents of matrices with positive entries, which can be effec- 
tively approximated [19|, though such permanents do not appear in the quantum networks. 



bosonic unitary networks, which are now under development |10l-ll2[|. 

It is beheved that the hardness of computing the permanent of an arbitrary matrix 
is related to the matrix rank. Different limits are possible in this respect, which can be 
roughly divided by the relation between the number of bosons A^ (i.e. the matrix size) and 
the number of modes M in the network (i.e. the maximum rank of appearing matrices). 



For instance, it is shown in Ref. J2l| that even an approximate computation of matrix 
permanents appearing in the bosonic linear networks is #P-hard in the limit oi N ^ oo 
and M ^ A^^ (this limit is related to the "boson birthday bound", i.e. the probability of two 
bosons landing into the same mode should be negligible). The computation of the matrix 
permanent of an arbitrary (complex) N x A^-dimensional matrix with the maximum matrix 
rank of -Rtv would apparently require on the order of A^^^ flop operations on a classical 



computer (see, for instance, Ref. |22[ |) 



On the other hand, no analytical or semi-analytical approach to approximate the matrix 
permanents giving the bosonic probability amplitudes was ever proposed. This is especially 
lacking in the limit of a large matrix size (i.e. A^ ^ 1), when the computational complexity 
makes the numerical evaluation practically hard. First of all, one would not expect the 
existence of an asymptotic analytical approximation in the limit A^ — ?> oo and M > N: 
it would then contradict the #P-hardness of numerical approximation of the permanent 
in this limit (since one could run a numerical scheme mimicking the analytical approach). 
Therefore, an asymptotic analytical approximation could exist only for M ^ A^. Such an 
asymptotic in 1/A^ analytical approach is developed below by reducing the summation of 
the Laplace expansion of the matrix permanent to a multidimensional integral of the saddle- 
point type in the limit A^ — ?■ oo and M being finite. The very existence of an asymptotic 
approximation is important and stems from the fact that it is not dependent on the matrix 
size N (i.e. the number of bosons in the network), but only on the solution of some 2M — 1 
bilinear equations giving the saddle points of the integral. 

Practical usage of the method depends on effective solution of two other problems of 
different type. One is the unitary matrix scaling problem, i.e. for a given unitary matrix U 
to find all diagonal scaling matrices X and Y with the complex-valued entries such that the 
product matrix XUY has given row and column sums (these matrices are the saddle points 
of the integral). The other problem is to derive the general formula for multidimensional 
saddle point method, i.e. in the case when the saddle points coalesce. Only the case of 



simple saddle points has a general solution so far [23| and a special case of two coalescing 



saddle points on the real axis is treated. [2j] However, the coalescing saddle points are 



exceptional points rather than the general case. In the case of simple saddle points, which is 
the general case, an explicit formula for the bosonic amplitudes is derived below. As we use 
the well-known Stirling formula for the factorial in one of the steps in the approximation, the 
approach is restricted to the case of inputs and outputs with no mode without at least one 
boson in it (though this latter limitation can be lifted, for simplicity sake it is not pursued 
in this work). The formula for the bosonic amplitudes is tested on the case of two-mode 
network (i.e. the beam splitter) where it shows a very good accuracy within the above 
described limitations. 

The rest of the text is organized as follows. In section [TTl an integral approximation for 
the bosonic probability amplitudes in the unitary linear networks is derived. The integral 
is evaluated by the saddle point method and an explicit formula is derived when the saddle 
points are simple. Some details of the calculations are relegated to [B] and O For complete- 
ness, a short derivation of the bosonic amplitudes as the matrix permanents is given in El 
In section IIIII the general formula is tested on the simplest case of the two-mode network 
(the beam splitter) where the sources of the error are identified and discussed. Section IIVI 
contains summary of the results and perspectives. 

II. THE SADDLE POINT METHOD FOR PROBABILITY AMPLITUDES IN 
LINEAR BOSONIC NETWORKS 

We consider the transition probability amplitudes of A^ bosons between the M input 
(l/i), . . . , I/a/)) and M output modes (I qi), . . . , Is'a/)) of an unitary linear network, which 



are given by (see, for instance, Refs. |20l. l25|) 



I , per(f/[ni,...,nAf|mi,...,mA/]) ,^. 



where the matrix f/[ni, . . . , nM\fnii • • • , fnu] consists of repeated rows and columns of the 
network matrix U (the order being insignificant), with total of Uk of the fcth row and total 
of mi of the /th column of f/, satisfying Ylik=i ^k = J2k=i "^^ ~ ^ (^^'^ more details, seeRll. 
By a successive application of the Laplace expansion for permanents [ij] one can express 
the boson probability amplitude ([1]) as an average over the lattice of M x M-dimensional 



matrices Ski {Ski ^ {0, 1,2,3.. .}) with given row and column sums, equal here to the 
distribution of bosons in the input and output modes, respectively: "^i^i Ski = nk and 
'^k=i ^ki = ^TT-i- Indeed, the first application of the Laplace expansion consists of dividing 
the matrix U[ni, . . . ,nM\iTii, . . . ,mM] into two parts, and the permanent into a sum of 
products of two permanents, one involving Ui first rows (all equal) and the other involving 
the rest of the matrix rows. The summation runs over all possible distributions of the 
column indices between the two factors in each product of permanents, such that one part 
contains rii column indices {wi, . . . , Wm, with ^n of them being equal to 1, S12 being equal 
to 2, etc) and the rest oi N — rii column indices {wn-^+i, • • • , wjy) go to the other part. Thus, 
application of the Laplace expansion as above described gives 



per([/[ni,...,nM|mi,...,mM]) = ^ per{U[ni\Su, 

1<wi<...<wn<N 

xper(f/[n2, ..., 71^1^1 - Su, ..., tum - Suj]) 



Sim] ) 



mil 



i^ii 



s s ,JrSu\{mi-Su). 

011,...,0iM t — 1 

xper(f/[n2, ...,nM\mi- Sn, ..., uim - Sim]), 



(2) 



where we have used that permanent of the first submatrix is equal to nil Y[i=i Un^' and, due 
to the permutational invariance of the matrix permanent, the r.h.s. contains functions of the 
total numbers of the repeated columns and not of the column indices themselves. Repeating 
the above procedure, by taking out successively each set of the repeated rows, we obtain^ 



peT{U[ni,...,nM\mi 

M 



.k=l 



Ski>0 



Skl^lT-l. 



.k=l 



,mM]) = 

M 



M 

]^nfc!mfc! 



.1=1 



^ jjSu 
^kl 



n 

k,l=l 



(3) 



Eq- © gives the permanent of a matrix having repeated rows and columns and has an 
interesting statistical interpretation.^ Indeed, the ratio of factorials which appears on the 



^ Eq. ^ could also be deduced from Eq. ([IJ and the formula for the bosonic probability amplitude derived 
in Ref. [261] by application of the Wick theorem; note that using the matrix permanent and the Laplace 
expansion is an impressive shortcut to the otherwise quite involved derivation. 

^ Besides a physical interpretation for the bosonic probability amplitude of Eq. ^■. the summation of the 
quantum probability amplitudes over all possible transitions through the linear network, i.e. a variant of 
R. Feynman's path integral formula. 



r.h.s. in Eq. ([3l) , divided also by A^!, is know as the Fisher- Yates distribution (see, for 
example, Ref. 27|]). It appears in Fisher's exact test of the applied mathematical statis- 
tics which uses the so-called contingency tables. The Fisher- Yates distribution gives the 
conditional probability of getting the actual matrix Sm (the contingency table) of the joint 
frequencies of two independent properties given the row and column sums (the margins) 
equal to their marginal frequencies (i.e. ni, . . . , um and mi, ... , tjim, respectively). It reads 
(using a shortcut notation {fj} for the set of indexed variables: fi,f2, . . .) 



'i-rAf ," 




"t-tM ," 



V{{S,,}\{n,,m,}) = ^ i^ \ (4) 

Note that the delta symbols in Eq. (J3]) restrict the summation to the contingency tables 
with given margins and precisely under these constraints the probabilities of Eq. (jl]) sum 
to 1. Then the matrix permanent of Eq. (JSj) is nothing but multiplied by A^! the value 
of the characteristic function x{{^ki}\{nk,fni}) = {^w{J2k i=i ^kiSki}) of the Fisher- Yates 
distribution at the parameters equal to the logarithms of the entries of U: 




per{U[ni, . . . , nM\mi, . . . , tum]) = N\ 

= N\x{{HUM)}\{nk,mi}). (5) 

Eq. (|5]), however, cannot be used for the numerical evaluation of the matrix permanents 
since the number of the contingency tables scales exponentially with A^ (nar nely , it scales 
exponentially with the margins for a fixed table size, see, for instance, Ref. [28|). On the 
other hand it is an indication to a possibility of an asymptotic approach: since the lattice 
of contingency tables is exponential in A^, if we divide each of these matrices by A^, then, as 
A^ — )■ oo, the resulting matrix p, pki = Ski/N, will belong to a dense lattice in the convex 
set of matrices with real-valued entries and given row and column sums. Therefore, the 
summation over such a lattice can be replaced by a multidimensional integration with the 
integration variables being p^i- Moreover, a large parameter, A^, appears in the exponent of 
the integrand allowing for an asymptotic evaluation of the matrix permanent. This is the 
approach pursued in the following. 

First of all, we need to approximate the Fisher- Yates probability (jl]) by a manageable 
smooth function of the integration variables pki- Using the approximate formula for the 



multinomial coefficient, given by Eq. (1B3P of [Bl we obtain for nkjiTik, Ski >!'■ 

1 

nM rik rrik 
k=l N N 



V{{Ski]\{nk.mi}) = {2txN)-'-^^ 



{M 



nM 
kI=lPkl 



xexp{-iVX({pfc,})}(l + 0(iV-^) 



where we have denoted by X the mutual information function, namely 



X{{Ski/N}) = n{{nk/N]) + n{{mk/N})-n{{Sui/N}) 



(6) 



M , 

^ Pki In I 



Pki 



kl=l 



\{nk/N){mi/N) 



(7) 



with the Shannon entropy function denoted by l-L. 

The summation in Eq. ([3]) can be replaced by integration as A^ — )■ oo, since the difference 
between the entries of two neighboring matrices (the A^~^-scaled contingency tables) pki is 
of order 1/A^. In this way, the Kronecker delta symbols must be replaced by the A^~^-scaled 
Dirac delta functions: 

1 



Ap, 



kl 



N 



-> dp, 



'kl, 



S-r^M c „ — )■ 









(8) 



There are only 2M — 1 independent constrains in the product of 2M Kronecker deltas in Eq. 
(|3]), since the contingency table elements Ski sum to A^ and this value trivially gives both the 
sum of row sums and column sums. Hence, the integration domain is (M — 1) ^-dimensional. 
Using these observations and some elementary algebra we obtain from Eqs. ([3]), (E]), and 



(M-l) 

per(?7[ni, . . . , n^lmi, . . . , niM]) ^ A^! -- 

ZTT 



-TT nkTIlk 



X / dfi{{pki})- 



exp < ~N 



M 

\| fc=i 

AiPki}) - Y.k,i=iVki In Uki\ I 



iM 



nM 
kl=lPkl 



(9) 



Here we have introduced the integration measure dfi{{pki}) over the (M — l)^-dimensional 
subspace in the convex set of all matrices p constrained by the total of 2M — 1 Dirac delta 
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functions from Eq. ([8]): 



.4 



dfi{{Pki}) 



M 



X 



.k,l=l 
'M-1 



kl 



M 



M 



HMEp- 



.fc=l 



.1=1 



M 



(10) 



1=1 \k=l 

The approximation error committed by Eq. iQ is estimated to have the multiphcative 
order ~ 1/A^, since this is the order of our approximation of the Fisher- Yates distribution 
by a smooth function in Eq. ([6]), whereas approximating the finite sum by an integral brings 
also the error on the order of the difference between the two lattice points, i.e. Apki ~ 1/A^. 
The multidimensional integral in Eq. (Q is of the standard form used in the asymptotic 
expansion in powers of 1/A^ by the saddle point (a.k.a. the steepest descent) method. 
However, since the error of approximation of the finite sum by the integral in Eq. (Q is of 
order 1/A^, only the leading term of the asymptotic expansion is meaningful. 

The saddle points (matrices) p are found as the extremals of the multivariate function 
in the exponent of the integrand (see, for instance, Ref. [23l]). In our case the function is 
(f) = ^{{pki}) — J2k 1=1 Pki ^^Uki with, however, only (M — 1)^ independent variables out of the 
total M^ matrix elements pki- Using the Lagrange multipliers A^ and fj,i one can equivalently 
minimize the function 

M M 

J" = Ii{pki}) - ^ Pki In Um - ^ (Ayt + Hi) Pki- (11) 

k,l=l k,l=l 

Equating the differential of J-" to zero we obtain that the saddle points have the following 
general form 

Pki = XkUkiVi, (12) 

where the complex parameters Xk and yi are determined by satisfying the margins imposed 
on p by the Dirac delta functions in the measure ( ITOl) . In other words, the diagonal matrices 
X = diag(xi, . . . , xm) and Y = diag(?/i, . . . , i/m) solve the matrix scaling problem for the 
unitary matrix U: the scaled matrix XUY must have the row and column sums specified 



Note that an arbitrary subset of 2AI — 1 delta functions can be used, see alsolCl 



by the boson distributions in the input and output modes, i.e. 



M M 

1=1 k=l 



^XkUkiyi = -j^, ^XkUkiVi = -jjr. (13) 



At this stage, let us recall the general formula for the leading term of an n-dimensional 



integral, as given by the saddle point approximation, when the saddle points are simple |23| : 



where the summation is over the contributing saddle points Zj. The determinant in the 
denominator of Eq. flT^ is of the Hessian matrix, i.e. the matrix composed of the second- 
order derivatives of (f){z), taken at the respective saddle point. The saddles Zj are found 
by a deformation, as allowed by the analyticity of 0(-2), of the integration domain in the 
complex-valued extension space, such that the deformed domain is contained in the steepest 
descent regions of the integrand. The next term in the asymptotic expansion, as compared 
to the leading term of Eq. flT^ . has the relative order of 1/A^. 
In our case the function reads 

M 

(j) = Z{{pki}) - ^ Pki In Uki (15) 

k,l=l 

and the matrix of the second derivatives with respect to the free (M — 1)^ variables from 
the M^ elements of the matrix p must be used. Therefore, an extension of the saddle point 
method to the constrained integration is needed, which runs as follows. We rewrite the 
constrains on the variables pki in Eq. flTUl) as a set of linear equations by introducing the 
matrix Cj^ki, where the enumeration order of the double index {k, I) is as follows {k, I) = 
{(1, 1), . . . , (1, M), (2, 1), . . . (2, M), . . . , (M, 1), . . . , (M, M)}, i.e. the index k runs slower 
then the index /. Then the constrains are expressed as follows^ 



7 ^ CjMPkl — Cj , Cj = { 



rii um mi niM-i 



ik,i) 

Cjki = < ' ~ (16) 

5j.M,i, M + l<j<2M-l. 



The specific subset of 2M — 1 constraints is in accord witfi tlie seiected measure in Eq. (ITU)) . 
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The matrix C in Eq. ( !T6|) has rank equal to 2M — 1. It can be partitioned into a nonsingular 
submatrix C*-^-* of the size (2M — 1) x (2M — 1) and the residue submatrix C^^^\ These 
induce a similar partition of the elements pki, treated as a vector with double index {k,l). 
Using the vector notation p, we can rewrite the system of constraints given by Eq. fITBl) as 
follows 

First, by introducing two vectors, one consisting of 2M — 1 dependent integration variables, 
^, and the other of (M — 1)^ independent ones, f], as follows ^ = C^^'p^^' and fj = p^^^\ we 
satisfy the constraints by fixing the value of S, according to Eq. (ITTl) (i.e. by integrating 
over ^ using the Dirac delta functions in Eq. flTU]) ) and obtain the rest of the measure d^ as 
follows 



n ^"^i 



(18) 



Next, due to the linearity of the constrains (1171) . the determinant of the matrix of second 
derivatives of (t>{{pki}) with respect to the independent variables can be evaluated from the 
full matrix of second derivatives with respect to all variables pki and the constraint equation 
(ITTl) . We get the following result 



det 



(- 



&" 



det<'[5,/]('S 



B 
I 



\d{p^^^^f, 
B = - [C^^))-' C^''\ (19) 

Here we have used the block matrix [B, I] constructed from the transposed (2M — 1) x 
(M — 1) ^-dimensional matrix B and the (2M — 1) x (2M — l)-dimensional matrix unit /. 
Furthermore, the determinant on the r.h.s. of Eq. ( TT9l) can be further simplified by using an 
identity which generalizes the well-known Sylvester's identity for the determinants of block 
matrices (see O for details) and is valid for a nonsingular matrix A: 

B 



det(C(^)) det<^ [BJ]A 



det(A)det CA-^C , 



(20) 



where the matrix B is as in Eq. (1191) . In our case the Hessian matrix A appearing in the 
second differential of = X{{pki}) — Ylik i=iPki ^^ Um is diagonal, i.e. 

M 



d^<P{{Pki}) =Y] — (dpki) 



(21) 



k,l=l 
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thus Eq. (120!) applies when there is no saddle point located at zero, i.e. for pki 7^ 0. Note 
that the exponent with the large parameter N in the integral on the r.h.s. of Eq. (Q taken 
at the saddle point p such that some of pki = would be infinite (thus this case is ruled 
out). Indeed, using the constraints flT^ . we get at a saddle point p^i = XkUkiHi'- 

exp j-iV {l{{pth) - E Pti ln/7,, ) \ = \[ 



nk \ mk^ 



k,l=l / J fe=l ^ '^ 

Now, using Eqs. dH]), and ([IB])- ([22]) into Eq. © and noticing that 

''t) ^ TT ^. 

Pkl 



Nvk 



(22) 



det 



dp^ 



k,l=l 



we obtain the final formula for the leading term approximation to the matrix permanent in 
the case of simple saddle points (our main result) 



peT{U[ni, . . . ,nM\mi, . . . ,mM]) ~ N\, 



M 



M 



nkfUk 



nuk iibk sr^ 

\ fc=l s 



nti 






Nx 



.M 



Nvi'^ 



Vdet(D'(p(«))) 



(23) 



» _ J')i 



,(«) 



where the summation over all contributing saddle points pl/ = x], UkiVi i^ implied and the 
matrix D' in the denominator is as follows 

^2 -\ -1 



D' = C 



d^ 



^7^2 



c. 



(24) 



with the matrix C given by Eq. ( 1T6|) . The sparsity of C allows one to easily find the explicit 

form of the (2M — 1) x (2M — l)-dimensional matrix D' . The simplifying observation is that 

the resulting determinant det(i5'), appearing in the denominator on the r.h.s. of Eq. ( 123|) . is 

actually equal to any of (2M — 1) x (2M — l)-dimensional principal minors (i.e. obtained by 

striking out the same column and row) of the full 2M x 2M-dimensional matrix D, defined 

as follows 

\ 



/ Hi 

' N 



D 



N 



p 



p 



mi 

N 



N 



(25) 
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(in Eq. (125|) we have used the guiding sohd hues to emphasize the block structure of D). 
Any of the above specified principal minors can be used due to the fact that any subset 
of 2M — 1 constraints from the full set of 2M ones could be used in the definition of the 
integration measure dfi of Eq. fITU]) . The determinant itself appearing in Eq. fl2^ can be 
further reduced to a much simpler form (see O, where the above described property of the 
principal minors of D is also directly verified). We have 

M 



det{D') 



.k=l 



■ M 

n 



rrik 



N 

k=l 



det(A'2-p'ArV) 
det(A;-p'A2V), (26) 



where we have denoted Ai = diag(^, . . . , ^), A2 = diag(^, . . . , ^), and the matrices A'^ 2 
and p' appear in D' (i.e. one of them has only M — 1 elements on the diagonal). Note that all 
the above discussion implies that the exact symmetries of the bosonic probability amplitudes 
are preserved in the approximation given by Eqs. ([1]), f fT2|) . fl23l) . and f l25|) . For instance, 
the inversion symmetry: (m^^, . . . , mgj^,j\nf^, ..., n/^,) = (nj,, . . . , n/^Jm^,, . . . , mg^^)*.^ 

III. EVIDENCE OF A GOOD ACCURACY OF THE SADDLE POINT APPROX- 
IMATION FOR THE BOSONIC PROBABILITY AMPLITUDES 

To apply the approximation ( 123|) we first have to solve the matrix scaling problem f lT3|) 
which is a bilinear system of equations in x and y. One can reduce the number of variables 
by half resolving one of the constraints in Eq. flT^ . for instance, for y. We obtain yi = 
J2k=i Ukii^k/Nxk) and the reduced system reads 

M 

j:-^Uk,U;,^ = ^^ / = 1...,M-1, (27) 

k,q=l " 

where there are M — 1 independent variables: r^ = Xk+i/xk-, /c = 1, . . . M — 1. The system 
(!27|) can be further rearranged into a polynomial in r^ system of equations with the highest 
order equal to 2(M — 1). It proves to be hard to solve analytically, even in special cases, such 
as of the Bell multiports (despite considerable efforts, the general solution was not found). 



^ This is also manifested by the exact canceUation of the probabiUty amphtudes, i.e. the generaUzed HOM 
effect, in Fig. 1 of section IIIII below. 
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To test the accuracy of the approximation (123|) we thus consider the simplest case of the 
beam-sphtter, i.e. just two modes M = 2. In this case there is one independent variable 
and Eq. (127|) is a trivial quadratic. 



As is known, 29|] the two and three dimensional networks are uniquely defined by the 
moduli of the unitary matrix elements \Uki\, whereas the 2M — 1 phases are scaled out by 
changing the unimportant phases of the basis states giving the input and output modes. 
Let us consider the sjTumetric beam splitter, which is given by the following matrix (as the 
phases are unimportant) 

t/ = I v^ ^ I . (28) 

\ 75 75/ 
Then Eq. (I27|) leads to the following quadratic equation 

^(^Y -2^,1^^ + 1 = 0, ^ = ^^-^ , (29) 

ni \X2/ y n2X2 2y/nin2 



7 




Thus, there are two saddle points p whose x components in Eq. flT2|) are 

e'-^ = 7 ^ zv^l-T^. (30) 

The y component is given by y, = (-^^ + ^^)/V2 and y^ = {^^+ ^J/V2. In Eq. m 
we have introduced the phase 0, however, it is a real value only under the condition that 
7^ < 1, i.e. when 

(An)" + (Am)2 < A^^, (31) 

with An = n2 — ni and Am = m2 — mi. Under the condition ( 13T1) the y component can be 
given as follows 

yi = sj^e^''^^y\ y2 = ^e^^'-'^y^ 



2y/m^m^' 



iS 



I 



AnAm ± ^JN"^ - {Anf - (Am) 



y/n^n^m^m^ 



(32) 



^ To obtain xi_2 from r^ = a;2/a;i we have taken into account that a multiphcation of the vector a; by a 
complex number does not alter the solution ri and, hence, it corresponds to the same saddle point p. 
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Since 4mim2(l — cr"^) = 4?2in2(l —7^) = A^^ — (An)^ — (Am)'^, the same threshold condition 
fl3T]) apphes to the phases ip and 6. When the condition fl3T|) is violated, the phases (p, ip, 
and 6 become complex valued (0 and ip become imaginary). This corresponds to a phase 
transition in the behavior of the corresponding probability amplitude of the beam-splitter 
(see below). 

The elements of the matrix p, i.e. the saddles of the integral in Eq. (123|) . can be obtained 
from the vectors x and y by Eq. flT2l) . We get: 



Pll 


i 

~ 2 






1 


^nl 


P12 


~ 2 


u 


P21 


1 
~ 2 


fr22 

[n 



+ 



+ 



^2 

^-=2liV + 



nin2 

iWV 

ni?22 
iWV L 



-7±zVl 



7^ 



7^ 



ni?22 



-7=Fivl -7^ 
7 ± iy 1 — 7^ 



(33) 
(34) 
(35) 

NNV'-' '■" (^^) 

which are valid for the whole square domain: | A?2| < A^, | Am| < A^. When the condition (13T|) 
is satisfied, i.e. when I7I < 1, both saddle points are complex valued and contribute to the 
integral in Eq. (!23|) . When it is violated, the two saddle points become real valued. However, 
one of them does not contribute to the integral, since it escapes out of the integration domain 
< Pfez < 1. Which one of the saddle points contributes depends on the sign of 7, i.e. the 
sign of Am, and the ratio n\jn2 (see also Fig. 1(c) below). 

Finally, after a simple algebra, the determinant given by Eq. fl26l) becomes 



det(D') 



T Tre" 



1- 


'An' 


2\ 2 


( 
1- 


'Am' 


' An 


2 


'Am' 


2\ 2 




_N _ 




. N _ 


/ 





(37) 



where we have used that 16nin2m,im,2 = [N"^ — (Ara)^) [N"^ — (Am)"^). 

Substituting Eqs. (!30l) . (!32l) . and ( !37l) into Eq. (!23l) and the resulting approximation into 
Eq. (II]) we obtain the saddle-point approximation for the bosonic probability amplitude of 
the symmetric beam-splitter ( l28l) . To compare with the exact result, we use the analytical 
expression for the bosonic probability amplitude for the case of U in Eq. (!28|) (see Ref. (3|): 



{mg,,mg^\nf^,nf^) = i^l^ 

q=0 



-1)'^ y/ni\n2\mi\m2\ 



q\{ni - q)\T{mi - q + l)r(m2 + q - ni + 1] 



(38) 
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where x is an overall phase. Eq. (138|) is convenient if the numerical software (e.g. the 
MATLAB) accepts negative integers as arguments of the F-function (using the F-functions 
is a simple way to incarnate the constrained summation appearing in the formula for the 
amplitude, since for negative integer n we have F(n) = oo). 

The comparison of the exact result 0381) and the saddle point approximation can be 
divided into three regions. The first region corresponds to the condition (13TD being satisfied. 
This region contains the generalized HOM effect for larger number of particles. [2, 13|] In 
this case the amplitudes of the two contributing terms in Eq. (123!) . corresponding to the 
two saddle points, have the same moduli and the only difference lies in their relative phase. 
The corresponding domain in of the two-dimensional plane of An and A?ti is the inside of 
the circle: [AnY + (Atti)^ ^ A^^. There is cancellation of the probability amplitudes for 

nn 

Til = n2 = N/2 and the odd values of mi. Fig. 1(a) (see also Refs. [2, ISj). In the saddle point 
approach, the cancellation is due to the symmetry of the two saddle-point contributions so 
that the only difference is their relative phase given by (—1)'"^ (note: the cancellation is 
captured exactly by the saddle-point approximation since it preserves the symmetries of the 
exact bosonic probability amplitudes). 

On the other hand, there is another regime: the exponential decay of the probability 
amplitude as mi approaches either or A^, see Fig. 1(c). This regime has not been reported 
previously (for instance, the approximation of Ref . [3|] only captures the oscillating regime). 
It appears when the condition (|3T1) is violated. In this case, there is just one contributing 
saddle point, the one which has the smallest moduli contribution to the permanent (this is 
similar to what occurs in the saddle-point approximation to the Airy function). 

The third region is about the circle (An)^ + (Am)^ ~ A^^. This region contains the 
coalescing saddle points and cannot be approximated by Eq. fl23l) valid for the simple saddle 
points only (their contributions diverge on this circle, which is due to the determinant (!37|) 
approaching zero). Outside this region, which is restricted to narrow neighborhoods of the 
points iTii = 10 and mi = 50 in Fig. 1(c), the saddle point approximation has a very 
good accuracy as is shown in Fig. 1(b), where the accuracy of the approximation ( l23l) is 
compared to that for the binomial coefficient, given by Eq. fIBSI) of [B] and used to build 
the approximation of the Fisher- Yates distribution ([6]). It is seen that the relative error is 
approximately twice as that of the approximation for the binomial coefficient. 

Finally we note that the transition from the two contributing saddle points, with the 
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0.15F: 



0.125- 




^ 0.03 
^ 0.02 



FIG. 1: Comparison of the saddle point approximation (shown by stars) with the exact result 
(shown by circles) for the bosonic probability amplitude of the beam splitter (I28p (to guide the 
eye, the numerical points are connected by lines). Panel (a): A^ = 30 and ni = 15. Panel (b) gives 
the relative error of panel (a) at the even points of mi (the squares) compared with the relative 
error of the used approximation of the binomial coefficient ( ) based on the Stirling formula (the 
circles). Panel (c): N = 60 and ni = 10, the two regions close to mi = 10 and mi = 50 are the 
vicinity of the circle (An)^ + (Am)^ = N'^ where the simple saddle point approximation of Eq. 
(j23p fails (diverges). For 10 < mi < 50 the two saddle points contribute, while for mi < 10 or 
mi > 50 just one saddle point contributes. 
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oscillating contribution (the probability amplitude as function of rrii in our case), to a single 
contributing saddle point, with an exponentially decaying contribution, is similar to the 
Airy function behavior, thus the related integral in Eq. f l23p can be, in principle, expressed 



m 



ts 



through a linear combination of the Airy function and its first derivative. Such resu. 
are available for the real-valued coalescing saddle points (for instance, in Refs. '. 
However, the method needs to be generalized to the complex-valued saddles before it could 
be applied to the integrals approximating the unitary matrix permanents. 

IV. CONSLUSION 

We have demonstrated that an analytical approach is possible for approximate evaluation 
of the matrix permanents giving the bosonic probability amplitudes in the unitary linear 
quantum networks. The approach is based on the observation that in the limit of a large 
number of bosons, N, in an unitary network of M modes with M <^ N, the matrix perma- 
nents giving the transition probability amplitudes of the bosons from the input modes to 
the output modes are intimately related to the characteristic function of the Fisher- Yates 
distribution giving the probabilities of the contingency tables appearing in the exact Fisher 
test in applied mathematical statistics. The large number of the contingency tables with 
large margins (i.e. for large A^) allows one to replace the finite sum over the dense lattice of 
the 1/A^-scaled contingency tables by a multidimensional integral having a large parameter 
A^ in the exponent of the integrand. The integral is then evaluated by an extension of the 
saddle point method to the constrained multidimensional integrals, developed in section [III 
In the case of the simple saddle points, an explicit formula for the matrix permanent is 
derived. 

The practical application of the method is conditioned on the solution of a matrix scaling 
problem, which reduces to a system of bilinear equations in 2M — 1 variables. The general 
analytic solution of the matrix scaling problem has not yet been found even for the simplest 
case of the symmetric multiport network, despite a considerable effort. However, the scaling 
problem has the A^-independent complexity in contrast to direct evaluation of the matrix 
permanents giving the bosonic probability amplitudes. 

The asymptotic approximation has been compared with the exact result on the example of 
the two-mode network, i.e. the beam-splitter, and is found to have good accuracy correlated 
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with the accuracy of the approximation for the multinomial coefficient used as a building 
block of the saddle-point approximation for the matrix permanent. Moreover, it was found 
that there could be various regimes of behavior of the bosonic probability amplitudes in 
unitary networks, which are dependent on the number of contributing saddle points. For 
instance, in the simplest case of beam-splitter considered here, there are two regimes: (i) the 
oscillating regime, when two saddle point contribute to the approximating integral and the 
generalized Hong-Ou-Mandel effects take place, and (ii) the regime of exponential decay of 
the probability amplitudes, when only one saddle point contributes. More complex regimes 
could be expected in the multimode unitary networks. 
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Appendix A: Bosonic amplitudes expressed as the matrix permanents 

Consider a linear network unitarily transforming A^ bosons from the input modes 
l/i), . . . , I/m) to the output modes \gi)^ ■ ■ ■ , Iqm), which are two orthogonal bases of the M- 
dimensional single-particle Hilbert space H. The quantum network is given by the M x M 
unitary matrix U, U'^U = I, which transforms the input modes into the output modes, i.e. 

M 

\fk) = Y.UM\9i). (Al) 

1=1 

The goal is to express the bosonic amplitude between the two Fock states |n/^, . . . , Ufj^^) and 
|mc,j, . . . , TUgj^), giving the input and the output states of N bosons: 

(A2) 



(A3) 



^M V^M 



|n/,,.. 


■^rrigj = , 




■ ? /jjv / ? 


rUg,,.. 


/ A! 


■■^93m)^ 



where X]fc=i ^fc ~ X]fc=i "^^ ~ ^^ ^^^ indices are defined 



as 



(^1, . . . ,^;v) = (W^' 2^_^, . . . ,M,.^., M). (A4) 

ni n2 riM 
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We have used on the r.h.s. of Eqs. (lA2p and (1A3|) the (unnormahzed) projection of the state 
of A^ particles in the tensor product of A^ single-particle Hilbert space if ® if (g) ... if on 
the symmetric (i.e. bosonic) subspace, given by the summation over all permutations r of 
the A^ indices of single-particle states divided by the number of all permutations, e.g. 

l/n, ■ ■ ■ JiN) = J^Yl I^Mi)) ® • • • ® l/vtiv,)- (A5) 

T 

The transition amplitude between the Fock states of Eqs. flA2p and (IA3I) then becomes 
the double summation over permutations of the indices in the / and g states in their inner 
products (i.e. the network unitary matrix elements). This is then converted to a single 
summation over all permutations of the column indices of the product of A^ unitary matrix 
elements Um by transferring one of the two permutations to the co-product indices, i.e. 

(M 
YlrnkW- 
k=l 

^J^_^^(9jAl)\firil)) ■ ■■■■ {9jAN)\fir(N)) 
<J T 

M V'" I 

k=l / cr T 

M \ ~2 

JJmfclnfc! ^ Ui,^j„,^^^ ■...■ ?7.^j^,(^), (A6) 

^fc=l / cr' 

where a' = a ■ t^^ runs over all permutations of A^ elements. One immediately recognizes 
the matrix permanent on the r.h.s. of Eq. (IA6p . where the matrix consist of repeated 
rows and columns of the network unitary matrix U, where the order is insignificant due to 
permutational invariance of the matrix permanent. Therefore, we introduce the notation 
U[ni, . . . , nAf |mi, . . . , rriM] for such a matrix and obtain the resulting amplitude proportional 
to the permanent of this matrix as in Eq. ([1]) of section [Til (see also Refs. 
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25|). 



Appendix B: Approximating the multinomial coefficient 

Here the possibilities for approximation of the multinomial coefficient, and hence the 
Fisher- Yates distribution, are detailed. We have the exact result for n > 1: [30 1 



n! = v/27r(n + ft„)(^)", (Bl) 
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where the free parameter ^„ is tightly bounded: 1/6 < 9n < 0.177. Interestingly, Eq. (jBTj) 
can be extended to all n > by carefully defining the limit 0° = 1 and redefining the lower 
limit on the free parameter as ^o = ^ < |- Then, after simple algebra, the multinomial 
coefficient becomes 

N\ _ exp {Nmrn)) ( 1 + f 

rM ■ ~ /77^ .rx»,f -I I _»^ r. All- \^^) 



Utinkl y^i2nN)M-i I j^ 



M 
k=l 



2k _|_ ""k 

N ^ N 



Here we have introduced the expected Shannon entropy function Hd^}) = 
— J2k=i 7^ 1^ ilv)- Since the free parameter 9 in Eq. (JB2I) is always divided by A^ ^ 1, the 
simplest approximation is obtained by dropping it out completely, thereby committing an 
error of order (9(A^~^) and restricting ourselves to n^ > 1. Assuming the latter, we obtain 

AT! exp(A^H({^})) 



nfli^fc! J(2vrAr)M-i nf=i t 



il + 0{N-^)). (B3) 



Finally, we could get an even better approximation for the multinomial coefficient (and 
for n > 0) by setting an uniform 6n in Eq. (IB2D . for instance 6n = 1/6, i.e. similar as 
in Gosper's approximation for the factorial. [SlJ Such a formula would then improve the 
approximation of the Fisher- Yates distribution by a smooth function in section [TTl However, 
for simplicity of the integral representation, we do not pursue this approach in the present 
work. 

Appendix C: Proofs of the determinant identities of section [III 

Let us first proof the generalization of Sylvester's determinant identity, i.e. Eq. ( 120|) of 
section [TTl To this goal, one can use the following auxiliary integral (where, for convenience, 
we replace Cj^ki by Cji, i.e. the double index by a single one, and denote m = 2M — 1 and 
n = M^) 

/m / n \ 

t/"xe-^^^'nW5^C',,a;, . (CI) 

Here x = (xi, . . . , a;„), det{A) ^ 0, and rank(C) = m. The determinant identity f[2U]) follows 
if we evaluate J by two different methods. We assume first that the real part A^i > in the 
Hermitian decomposition oi A: A = A^ + iAj. The first method is direct integration over 
the n — m independent variables y extracted from variables x. Suppose that the full rank 
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submatrix C^^^ is given by the first m columns of C and C'^^^^ is tlie remaining submatrix. 

Ttien y = {xm+i, ■ ■ ■ ,Xn)- Resolving tlie constraints given by tlie Dirac delta functions in 

Eq. (]C1|) we obtain the dependent variables x^^^ = By, where B = —(C^^^)^^C^^^\ and the 

whole vector 

^ B 

I 
Then the integral in Eq. fICip becomes 



X 



y- 



(C2) 



J 



1 



|det(C(^))| 



r-'-yexp{-y[B,I]A 



B 
I 



(C3) 



where the real part of the quadratic form in the exponent is positive definite. The Gaussian 
integral in Eq. fICSP can be easily evaluated and we obtain 



J = 7r'^|det(C(^))rMet [B,I]A 



B 
I 



(C4) 



On the other hand, one can also use the Fourier representation for the Dirac delta functions in 
the integrand of Eq. (1C1I) and, by interchanging the integration order, evaluate the integral 
J on the whole x space and then take the inverse Fourier transform. As all appearing 
integrals are Gaussian they are easily evaluated. We have 



J 



R^ 



(F'x I ^ exp{— xAx + iXCx} 

(27r)-yd^^ ^^ 4 ^ 



n — m 

TC 2 



det(A)det(CA-^C 



(C5) 



Comparing the two results in Eqs. (IC4p and (ICSp we obtain the determinant identity 
of section [Til valid for the nonsingular matrices A with the positive definite real part. The 
validity can be extended to arbitrary nonsingular matrices A by uniqueness of the analytic 
continuation, observing that the r.h.s.'s of Eqs. ( lC4p and flC5P are analytic functions of the 
entries of A. 

Now, let us verify that all the (2M — 1) x (2M — l)-dimensional principal minors of the 
matrix D fl2^ . i.e. obtained by striking out one row and one column with the same index. 
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have the same determinant. This is the consequence of the unique null eigenvector v of the 
matrix D, Dv = 0, in the form v = [I, — /], where / is the M x M-dimensional matrix unit 
(which fact follows from the definition (1251) and that rank(C) = 2M— 1). Consider now the 
adjoint matrix D composed of the minors of D (in general, the adjoint of A is the matrix A 
satisfying A A = det(A)/). The principal minors of D are equal to the diagonal entries of D, 
thus we have to verify that the diagonal elements of D are equal. We have: DD = which 
is possible only if there is a vector u such that D = uv. Since D is a symmetric matrix, such 
is also D and we have D = avv for some scalar a. The latter equality means Djj = a and 
hence the (2M — 1) x (2M — l)-dimensional principal minors of D are all equal. This fact 
allows us to use the root of any of the (2M — 1) x (2M — l)-dimensional principal minors of 
D in the denominator of the r.h.s. in Eq. fl25]) . 

Finally, let us obtain the value of the (2M — 1) x (2M — l)-dimensional principal minors 
oi D. To this goal the following determinant identity of the block matrices can be used 

det I ^ ^ I = det(Ai)det(A4 - A3A-M2) 
A3 A4 



det(A4)det(Ai - Aa^J^Ag 



(C6) 



which is a block-matrix generalization of the rule for determinant of 2 x 2-dimensional 
matrices. Extracting one of the (2M — 1) x (2M — l)-dimensional principal minors D' from 
Eq. fl2Sl) we obtain 

rik 



det{D') 



fe=i 

M 

n 



k=l 



N 



det(A'2-p'ArV) 
det(A;-p'A2y) 



(C7) 



where we have denoted Ai = diag(^, . . . , ^), A2 = diag(^, . . . , ^^), and the matrices 
A.'^ 2 ^^'i P' appear in D' . 
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